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^ ABSTRACT 

Quasi-Periodic Oscillations (QPOs) observed at the tail end of Soft Gamma Repeaters giant flares are com- 
monly interpreted as the torsional oscillations of magnetars. From a theoretical perspective, the oscillatory 
motion is influenced by the strong interaction between the shear modes of the crust and magnetohydro- 
^ dynamic Alfven-like modes in the core. We study the dynamics which arises through this interaction, and 

' ' present several new results: (1) We show that discrete edge modes frequently reside near the edges of the core 

Alfven continuum, and explain using simple models why these are generic and long-lived. (2) We compute 
the magnetar's oscillatory motion for realistic axisymmetric magnetic field configurations and core density 
["T^ profiles, but with a simplified model of the elastic crust. We show that one may generically get multiple 

gaps in the Alfven continuum. One obtains strong discrete gap modes if the crustal frequencies belong to the 
gaps; the resulting frequencies do not coincide with, but are in some cases close to the crustal frequencies. 
(3) We deal with the issue of tangled magnetic fields in the core by developing a phenomenological model 
to quantify the tangling. We show that field tangling enhances the role of the core discrete Alfven modes 
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O and reduces the role of the core Alfven continuum in the overall oscillatory dynamics of the magnetar. (4) 

We demonstrate that the system displays transient QPOs when parts of the spectrum of the core Alfven 
^ modes contain discrete modes which are densely and regularly spaced in frequency. The transient QPOs 

I— I are the strongest when they are located near the frequencies of the crustal modes. (5) We show that if the 

^ neutrons are coupled into the core Alfven motion, then the post-flare crustal motion is strongly damped and 

^ has a very weak amplitude. We thus argue that magnetar QPOs give evidence that the proton and neutron 

QQ components in the core are dynamically decoupled and that at least one of them is a quantum fluid. (6) We 

show that it is difficult to identify the high-frequency 625Hz QPO as being due to the physical oscillatory 
CO mode of the magnetar, if the latter's fiuid core consists of the standard proton-neutron-electron mixture and 

is magnetised to the same extent as the crust. 

^— «) Subject headings: Neutron stars 



1. Introduction netic energy powers the flares (Duncan 1998). This opens a 

unique possibility to perform an asteroseismological anal- 
. ^ Since the discovery of quasi periodic oscillations (QPOs) ygig neutron stars, and possibly obtain a new observa- 
^ in the lightcurves of giant flares from soft gamma repeaters ^-jo^al window to study the neutron-star interiors. Many 
in (SGR) (Israel et al., 2005; Strohmayer & Watts, 2005; authors have considered torsional modes to be confined 
^ Watts & Strohmayer, 2006; Barat et al., 1983) there has ^ magnetar crust, and have shown that seismological 
been considerable interest in their physical origin. One of information about such modes would strongly constrain 
the appealing explanations is that the QPOs are driven physics of the crust (Piro 2005, Watts & Strohmayer 
by torsional oscfllationsQ of the neutron stars whose mag- 2OO6, Watts & Reddy 2007, Samuelsson & Andersson 2007, 
Steiner & Watts 2009). 

^ By torsional oscillations we mean those which are nearly incompress- . . , , 1 ; 1 1 1 1 

ible. Modes with compression haye high restoring forces and feature However, it was quickly understood that the thcorctl- 

much higher frequencies than most of the obseryed QPOs. cal analysis of magnetar osciUations is complicated by the 



1 



presence of an ultra strong magnetic field (B ^ 10^'^ — 10^^ 
G) that is frozen into the neutron star and penetrates both 
the crust and the core. The field provides a channel for 
an intense hydro- magnetic interaction between the motion 
of the crust and the core, which becomes effective on the 
timescale of <C 1 second (Levin 2006). Since the QPOs are 
observed for hundreds of seconds after the flare, it is clear 
that the coupled motion of the crust and the core must be 
considered. In recent years, significant theoretical effort 
has gone into the study of this problem (e.g., Glampedakis 
et al. 2006, Levin 2007, Gruzinov 2008b, Lee 2008). This 
paper's analysis is based, in part, on an approach of Levin 
(2007, L07). 

To make progress in computing the coupled crust-core 
motion, L07 has studied the time evolution of an axisym- 
metric toroidal displacement of a star with axisymmetric 
poloidal magnetic field. In that case the Alfven-type mo- 
tions on different fiux surfaces decouple from each other, a 
well-known fact from previous MHD studies (for a review 
see Goedbloed & Poedts 2004, hereafter GP). One can then 
formulate the full dynamics of the system in terms of dis- 
crete modes of the crust which are coupled to a continuum 
of Alfven modes in the core. L07 demonstrated that (i) 
the global modes with frequencies inside the continuum 
are strongly damped via a phenomenon known in MHD 
as resonant absorption (see GP), and (ii) in many cases, 
asymptotically the system tends to oscillate with the fre- 
quencies close to the continuum edges. This result was 
later confirmed by Gruzinov 2008b, who has used a pow- 
erful analytical technique to solve the L07's normal-mode 
problem (Gruzinov noted that the resonant absorption is 
mathematically equivalent to Landau damping). Oscilla- 
tions near the continuum edge frequencies were also ob- 
served in a number of numerical gcnciral-rclativistic MHD 
simulations of purely fluid stars (Sotani et al. 2008, Co- 
laiuda et al. 2009, Cerda-Duran et al. 2009). 

Apart from finding QPOs near the continuum edges, 
L07's dynamical simulations identified transient QPOs 
with drifting frequencies; these were transiently amplified 
near the crustal frequencies. No explanation for the origin 
of the drifts was given. 

In this paper, we extend the previous analyses of the 
hydro-magnetic crust-core coupling in an essential way. In 
section 2, we re-analyse L07's toy model of a single oscilla- 
tor coupled to a continuum, and we show that this system 
generically contains the edge normal modes with frequen- 
cies near the continuum edges. We show that these modes 
dominate the late-time dynamics of the system, and de- 
velop a formalism which allows one to predict analytically 
the edge mode's amplitude from the initial data. We then 



explore the effect of viscosity on the system (introduced 
as a friction between the neighbouring continuum oscilla- 
tors), and show that the edge mode is longer lived than all 
other motions of the system. We also provide a non-trivial 
analytical formula for the time dependence of the overall 
energy dissipation. 

In section 3, we describe how transient QPOs, not asso- 
ciated with the normal modes of the system, are obtained 
when parts of the core spectrum consist of densely and 
regularly spaced discrcitc modes (and in section 5 we show 
that such an array of discrete modes is expected when the 
magnetic field in the core is not perfectly axisymmetric 
but has some degree of tangling) . As a by-product of our 
analysis, we explain the origin of the QPO frequency drifts 
seen in L07 simulations. We provide simple analytical fits 
to the drifts, and show that when the regularity of the con- 
tinuum sampling is removed (e.g, when the frequencies are 
sampled as random numbers picked from the continuum 
range), the drifts disappear. 

In section 4, we set up models with a more realistic 
hydro-magnetic structure of the neutron-star core. We 
show how to find the continuum modes and their coupling 
to the crust for an arbitrary axisymmetric poloidal field, 
with an arbitrary density profile on the core (L07s cal- 
culations, for simplicity and concreteness, were restricted 
to constant-density core and homogeneous magnetic field) . 
We treat a more complicated case of a mixed axisymmet- 
ric toroidal-poloidal field, with radial stratification, in the 
Appendix B. We demonstrate that for realistic field con- 
figurations, the Alfven continuum of modes coupled to the 
crust may show a number of gaps. If a crustal mode fre- 
quency belongs to one of these gaps, a strong global dis- 
crete mode arises which dominates the late-time dynamics 
and whose frequency also belongs to the gap. The fre- 
quency of the gap global mode docs not generally coincide 
with, but is often close to that of the crust. We suggest 
that it was these gap modes that appeared in Lee's (2008) 
calculations as well-defined discrete global modes. 

So far, only axisymmetric magnetic fields have been con- 
sidered in the magnetar-QPO literature, with the Alfven 
continuum modes occupying the flux surfaces of the field. 
In section 5 we argue that if the field is not axisymmetric 
but instead is highly tangled, then the Alfven continuum 
modes become localized within small regions of individual 
field lines, and therefore become dynamically unimportant. 
Instead, a set of discrete Alfven modes appears, with the 
spacing between the modes strongly dependent on the de- 
gree of field tangling. We devise a phenomenological pre- 
scription which allows us to parametrize the field tangling 
for computing the dynamically important modes, and in- 



2 



troduce an easily solvable "square box" model suitable for 
exploring the parameter range. 

Finally, in section 6, we use the suite of models built 
in the previous sections to explore their connection to the 
QPO phenomenology. We find that 

(a) within the standard magnetar model, it is possible to 
produce strong long-lived or transient QPOs with frequen- 
cies in the range of around 20 — 150Hz, but only if the 
neutrons are decoupled from the Alfven-like motion of the 
core; this implies that at least one of the baryonic compo- 
nents of the core is a quantum fiuid. 

(b) Our models could not produce the high-frequency 
625Hz QPO within the standard paradigm of a magnetar 
core composition. 

2. An oscillator coupled to a continuum: edge 
modes 

In this section, we study the motion of a harmonic os- 
cillator (which we hereafter call the large oscillator) which 
is coupled to a continuum of modesj^ This model was in- 
troduced in LOT and it provides a qualitative insight into 
the behaviour of crustal modes (represented by the large 
oscillator) coupled to a continuum of Alfven modes in the 
core of a magnetar. L07 found that if the large oscillator's 
proper frequency was within the range of the continuum 
frequencies, then the late-time behaviour of the system 
was dominated by oscillatory motion near the edges of the 
continuum interval. Here, we give an explanation of this 
phenomenon in terms of the edge modes. Our analysis al- 
lows us to use initial data and predict the displacement 
amplitudes and frequencies of the system at late times. 

The model consists of the large mechanical oscillator 
with mass M and proper frequency ujq, representing a 
crustal elastic shear mode. Attached to the large oscillator 
is a set of N smaller oscillators of mass m„ and proper 
frequency aj„ constituting a quasi-continuum of frequen- 
cies ujn (where n = 1,2,..., N). The continuum is achieved 
when iV — ^ oo while the total small-oscillator mass Sm„ 
remains finite. The convenient pictorial representation is 
through suspended pendulae, as shown in Fig. [l] (see also 
Fig. 2 of L07). 

The equations of motion arc obtained as follows. Each 
small oscillator is driven by the the motion of the large 



■^In many areas of physics similar models have been studied, notably 
in quantum optics and plasma physics. By contrast with the case 
studied here, in these models the range of the continuum frequencies 
is not limited. 




Large oscillator 



N small oscillators 



Fig. 1. — Schematic picture of the toy- model. A large number 
of small pendulae, representing the (quasi-) continuum, are 
coupled to one large pendulum, representing the crust. 



oscillator: 



(1) 



where a;„ is the displacement of the n'th small oscilla- 
tor in the frame of reference of the large oscillator, xg 
is the displacement of the large oscillator in the inertial 
frame of reference, and the right-hand side represents the 
non-inertial force acting on the small oscillator due to the 
acceleration of the large one. The large oscillator experi- 
ences the combined pull of the small ones: 



Mxo + MCoIxq = ^ 



■niiUJi X, 



(2) 



Here wq is the frequency of the big pendulum corrected 
for the mass loading by the small pendulae, i.e. Wq = 
Lol {M + m,) /M. 

2.1. Time-dependent behavior. 

In this subsection we explore the behavior of this sys- 
tem by direct numerical simulations. We found this to be 
helpful in the building of our intuition. We defer a semi- 
analytical normal-mode analysis to the next subsection. 
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Fig. 2. — Displacement of the big oscillator as a function of 
time. 

We follow L07 and for concreteness concentrate on a 
specific example; it will be clear that the conclusions we 
reach are general. We choose wq = Irad/second and mass 
M = 1. We choose a total number of 1000 small pendu- 
lae with frequencies ujn — (0.5 + n/1000)rad/second and 
masses to„ = m = 10"**, to mimic the continuum fre- 
quency range between O.Srad/second and 1.5rad/second. 
The simulation is initiated by displacing the large oscillator 
while keeping the small pendulae relaxed (this mimics the 
stresses in the crust), and then releasing. The subsequent 
motion of the system is then followed numerically by using 
a second order leapfrog integration scheme which conserves 
the energy with high precision. The resulting motion of the 
large pendulum can be decomposed into three stages (see 
Fig. [I and Fig. [3]): 

(1) During the first 50-60 seconds, there is a rapid expo- 
nential decay of the large oscillator's motion, during which 
most of the energy is transferred to the multitude (i.e., the 
'continuum') of small oscillators. This is the so-called phe- 
nomenon of "resonant absorption" , which has been stud- 
ied for decades in the MHD and plasma physics commu- 
nity (e.g., lonson 1978, HoUweg 1987, Goedbloed & Poedts 
2004, L07, Gruzinov 2008b). In this first stage, the am- 
plitude of the big pendulum motions drops by a factor of 
- 100. 

(2) After ~ 60 seconds, the exponential decay stops 
abruptly as the large oscillator now reacts to the collec- 
tive pull of the small ones. This second stage is character- 
ized by a slow algebraic decay of the amplitude of the big 
pendulum displacement. Gruzinov (2008b) explains this 
as being due to the branch cut in the oscillator's response 



Fig. 3. — A zoomed-in version of Fig. [2] The blue horizontal 
lines denote the theoretically predicted amplitude of the domi- 
nating upper edge-mode (see section 2.3). 

function. 

(3) The motion of the large oscillator stabilizes at a 
constant level (L07 missed this stage in his simulations, 
which he stopped too early). Fourier transform reveals 
2 QPOs at the frequencies close to the continuum edges, 
uj = 0.5 and lu = 1.5; the same QPO frequencies can be 
observed in the previous stage (2) as well. 

What is the origin of the QPOs, and how is this eventual 
stability established? In Fig. [4] and [5j we show how the 
amplitude of the small oscillators evolves with time. 

After the initial resonant absorption phase, the ampli- 
tude is distributed as a Lorentzian centered on the fre- 
quency around a; = 1; this is because the small oscillators 
in resonance with the large one are the ones which gain the 
most energy. However, in subsequent times we see that the 
energy exchange occurs between the small oscillator^ and 
that the net result of this exchange is the energy flow to- 
wards the oscillators whose frequencies are near the edges. 
By the time the third stage begins, the amplitudes of the 
oscillators near the edge stabilize and their phases become 
locked. They are pulling and pushing the large oscillator in 
unison. In the next subsection, we show that this behavior 
is due to the presence of the edge normal modes, and we 
shall derive their frequencies and amplitudes. 



This is much akin to the well-known phenomenon of resonant energy 
exchange between 2 equal-frequency pendulae hanging on the same 
supporting wall. 
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2.2. Finding eigenmodes 




0.6 1 

oscillator frequency 



Fig. 4. — The colored curves show the amphtudes of the smaU 
oscillators during the numerical simulation, at diilerent times 
t. 
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oscillator frequency 



Fig. 5. — A zoomed-in version of Fig. |4] At later times energy 
is transferred to the oscillators near the edge of the continuum. 



In this section we deal with the system of coupled har- 
monic oscillators, and one should be able to find its normal 
modes using the standard techniques (Landau and Lifshitz 
mechanics, §23). However, the fact that all small oscilla- 
tors are attached to the large one, and there is no direct 
coupling between the small oscillators, allows us a signif- 
icant shortcut (in Appendix A, we treat a more general 
problem of several large oscillators coupled to a multitude 
of the core modes). We proceed as follows: 

Suppose that we impose on the large oscillator a peri- 
odic motion with angular frequency £7, by driving it exter- 
nally with the force F^xt — FQ{n) exp{int) . This motion 
in turn drives the small oscillators according to Eq. ([T]): 



which has the steady state solution: 



(3) 



(4) 



where we have omitted the time dependent factor exp{int) 
on both sides. The combined force /cont of the small oscil- 
lators acting back on the large one (see Eq. ([2])) is given 

by 



/cont {^) ^^rrinUjl-^—^xo. 



According to Newton's second law, 

Fo m + /cont m = -Mi^f - ui)xQ. 



(5) 



(6) 



If f2 corresponds to the normal-mode frequency, then 
Fo{n) = 0. Hence by substituting Eq. ^ into Eq. ^ 
we get the following eigenvalue equation for 



G(r2) = M {ujI - 1]2) _ ^ 

n 

In the continuum limit N - 
becomes 



G{n) ^ M (wg - n^) - 



0. 



(7) 



cxD, the above equation 



j2 -02 



= 0, 



(8) 

where p(w) = dm/du) is the mass per unit frequency of the 
continuum modes. 

In the discrete case, the solutions of Eq. Q are — 1 
frequencies fi^ that are within the quasi-continuum (w^ < 
r^i+i < for i — l,2,...N — 1; 'quasi-continuum 
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modes'), and 2 modes with frequencies fliow and fihigh 
(the 'edge-modes') that are near the edges, but outside, 
of the continuum (i.e. fiiow < and fihigh > ^n)- It 
can be shown from Eq. ([7| that in the hmit TV ^ 1 and 
m„ <C M, the contribution of the small oscillator to the i- 
th quasi-continuum mode is completely dominated by the 
pendulae that are in close resonance with the mode. More 
precisely, one can show that as the number of oscillators 
A'^ increases and m„ decreases, the number of small oscil- 
lators contributing to the mode energy remains constant. 
However, for the two edge modes there is no such singular 
behavior in the limit of large N, and consequently they 
play a special role in the dynamics of the system. 

This last point is clearly seen in the continuum case 
represented by Eq. (|8|. The eigenvalue equation has no 
real solutions in the range of small-oscillator continuum 
Wmin < < Wmax, siucc the rcsponsc function G{n) is 
ill-defined in this intervaQ However, the edge modes on 
both sides of the continuum interval remain, and their fre- 
quencies can be found by numerically evaluating the zero 
points of G{fl) in Eq. ([8|. For the numerical of the pre- 
vious subsection, one finds fiiow — 0.5 — 8.2- 10^^ and 
f^high = 1.5 + 8.6- 10"^. Analytically, one can find the fol- 
lowing scaling for the distance SuJcdgc between the mode 
frequency and the nearest edge Wcdgc of the continuum 
range: 



SuJc 



C exp ■ 



M\nl 



edge I 



P(^cdgc) 



(9) 



edge 



where C is a constant of order unity. The larger is the 
density of continuum modes at the edge p{uJcdge), the fur- 
ther is the edge mode pushed away from the continuum 
range. It is particularly interesting to consider the case 
when the continuum interval is limited by a turning point 
(L07) with the divergent density of states near the edge, 
p{uj) = A/ ^/jti'^-a^dgolj where A is a constant. In this case 
the distance from the edge-mode frequency to the nearest 
continuum edge is given by 



Sui, 



edge 



Auj 



7/2 
edge 



Wedg 



M|^g-^?dgel 



(10) 



The quadratic dependence in Eq. ( 10 ) vs. the exponen- 
tial dependence in Eq. ([9| implies that the continua with 



^There is a complex solution if the integration in the expression for 
G(n) is performed along the contour chosen according to the Landau 
rule. One then obtains a "resonantly absorbed" or "Landau-damped" 
mode (Gruzinov 2008b, LOT), which exactly represents the exponen- 
tial decay of stage (1) in our numerical experiment of the previous 
subsection. 



turning points typically feature much more pronounced 
edge modes and stronger QPOs than the ones with lin- 
ear edges. In the next section, we show how to calculate 
the edge-mode amplitudes and QPO strengths from the 
initial data. 

2.3. Late time behavior of the system 

In the continuum limit, the only modes with real oscil- 
latory frequency are the edge modes. Thus, as we demon- 
strate explicitly below, they dominate the late-time dy- 
namics of the system when the number N of small oscilla- 
tors becomes large. Our analysis proceeds as follows: 
Lets define a new set of variables, expressed as a vector X 
with components Xq = y/Mxo and Xn — y/rn^(xo + Xn) 
for n — 1,...,N. With these new variables, the kinetic 
energy of the system is a trivial quadratic expression 



K =lx ■ X, 



(11) 



where the inner product of two vectors A and B is de- 
fined a.s A ■ B = YijLQAjBj. The potential energy is a 
positive-definite quadratic form, whose exact form is unim- 
portant here. The mutually orthogonal eigenmodes X^ can 
be found via a procedure outlined in the previous sectiorj^ 
Their eigenfreguencies 17^ are identified by finding zeros of 
G{Q) in Eq. ([t]), and the corresponding eigenvector com- 
ponents are given by 



X^ = 1 



(12) 



XL = 



Lets assume that we initiate our simulation by displac- 
ing the large oscillator by an amount xq (0) while keeping 
the small oscillators relaxed x„ (0) = and all initial ve- 
locities at zero. In the new variables, the initial state of the 
system is given by the vector X(0), where Xq = ^/MxoiO) 
and Xn = y/m^XQ^O). The time evolution of the system is 
given by: 

X{t) = T,n,cos{n,t) (x' ■ X'^ ^ (x{0) ■ X'^ X\ (13) 
Substituting the initial conditions, and the expression in 



Eq. ( 12 ) for the eigenvector components, we get 



X{t) 



Ecos(Qit) 
M 4- V — 



x\ 



(14) 



^Alternatively, they can be found by diagonalizing the potential- 
energy quadratic form. 



6 



The coordinate of the large oscillator is simply given by 
xoit)=Xoit)/VM. 

For the continuum of small modes, the above expansion 
breaks down, since the eigenvalue equation has no real so- 
lutions inside the continuum range. However, the edge 
modes are well defined, and they determine the dynamics 
at late times. Therefore, for the continuum case we can 
still write down the analogous expression which is valid 
only at late times: 



x{t) = Sfi^j^^ cos(rjcdgcO 



X{0) ■ Xedge 
edge * ^edge 



X. 



X, 



edge 



(15) 



The sums of Eq. (14 1 are replaced with the correspond- 



ing integrals, and we have the following expression for the 
displacement of the large oscillator at late times: 



E 



xait) = a;o(0) cos(f2cdgeO 



M + J du}p{uj) 



cdgc> 

(16) 

This expression is in excellent agreement with the numer- 
ical simulations. In the numerical example of subsection 
2.1, the upper edge mode dominates the late-time behavior 
of the system, and its calculated contribution is plotted in 
Fig. [3] together with the numerically simulated motion. 

2.4. The effect of viscosity 

We now add an extra degree of realism by introduc- 
ing viscous friction into the system. In MHD, continuum 
modes are spatially localized, and the effect of viscosity 
is to frictionally couple the neighboring modes (see, e.g., 
HoUweg 1987). In our simple model we introduce viscosity 
by adding frictional forces between the neighboring oscil- 
lators, 

/ri,ri+l = -/n+l,n = liXn " in+l), (17) 

where fn,n+i is the force from the n'th oscillator acting on 
the {n + iyth. 

We now calculate how the system dissipates energy as a 
function of time. We will show that it occurs in two stages 
(see Fig. [6|: (1) Initially, the small oscillators are strongly 
and simultaneously excited by the " Landau-damped" large 
oscillator, then they become dephased, with the average 
relative motion between the neighboring oscillators grow- 
ing linearly in time. This leads to a very rapid dissipa- 
tion of the bulk of the initial energy. (2) The edge modes 
persist, since the participating small oscillators move in 
phase and the energy dissipation is small. The energy of 
the modes is damped exponentially on a timescale much 
longer than that of the first stage. 



Viscous dissipation 




Fig. 6. — The red squares show the viscous dissipation of the 
total energy during the numerical simulation. The dotted blue 



curve shows the analytical solution from Eq. 1 23 1 



The dissipated energy is given by 



(18) 



In the continuum limit, the small oscillators are labeled 
not by a discrete index n, but by a continuous variable A. 
The expression for the dissipated energy is then 



Wdiss = / dXj 



dXdt 



(19) 



where 7 is the viscous coefficient. After the initial expo- 
nential damping of the large oscillator and the excitation 
of the small oscillators, the latter initially move indepen- 
dently, with 

xx{t) ^5:{\)cos[ujxt], (20) 

where x{X) is the amplitude of the A'th oscillator. From 
the above equation, we then obtain 



. dXdt 



1 



{[d{xxuJx)/d\f +u:lxl{dujxld\fe] 

(21) 

where the (...) stands for time-averaging over many oscilla- 
tion periods. For times t ^ dlogxx/ duj\ the second term 
on the right-hand side of Eq. (21 ) dominates. For a simple 

const, 



model with dujx/ d\ = const and p(w) 
dE/dt cx -At^E, 



(22) 



where E is the total energy of the system and A = 
/ p){dLOxl dX). The analytical solution for the energy and 
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O.B 1 1.2 

oscillator frequency 

Fig. 7. — As in Fig. Q, this figure shows the amphtudes of the 
small oscillators at different times t. The energy of most oscil- 
lators is drained due to viscous dissipation. At late times, only 
the oscillators near the edges of the continuum have substantial 
energy. 



dissipated power, 



E = 



(23) 



dE 
Itt 



= At^Soexp 



agrees very well with numerical simulations, see Fig. |6] 
While the equations above were derived for restrictive as- 
sumptions {d(jj\/d\ — const and p{ui) = const) , we found 



that the analytical formulae in Eq. (23) provide a good 
fit for a large variety of simulations. This is because it is 
the small oscillators with the frequencies near that of the 
large oscillator which carry most of the energy, and in that 
narrow band our approximations hold. 

After the energy dissipation due to dephasing is over, 
only the edge modes remain. This is illustrated in Fig. [?) 
where we show how the energies of the small oscillators 
evolve with time. At late times, only the oscillators taking 
part in the edge modes move substantially; this is because 
they remain in phase and do not dissipate much. At this 
stage the energy is drained by slow exponential decay of 
the edge modes. 

3. Transient and drifting QPOs 

Finite-size MHD systems feature a mix of continuum 
and discrete modes (see Poedts et al., 1985 and GP). For 



axisymmetric field configurations the continuum modes oc- 
cupy the whole flux surfaces and play an important role in 
the oscillatory dynamics; this was the motivation for L07 
and our study of the previous section. We will argue in sec- 
tion 5 that if the core field is highly tangled, the continuum 
modes become localized in space and discrete core modes 
will play a more important role. Thus it is important to 
study the case when the the crustal modes are coupled to 
a set of discrete core modes. In this section we show that if 
the frequencies of the discrete modes are regularly spaced 
in some frequency intervals, then the system displays tran- 
sient QPOs that are entirely missed by its normal-mode 
analysis. This is interesting from the observational point 
of view, since many of the observed magnetar QPO fea- 
tures are transient. 

Suppose that a set of discrete modes are located in the 
interval Ao; around frequency and are separated by a 
regular frequency interval 5u), and assume the following 
hierarchy: 

(5w < Aw < wq. (24) 

After the modes are excited, they are initially in phase but 
will de-phase rapidly on the timescale l/Aw. However, at 
times tn = 2Trn/Suj the modes come into phase again and 
pull coherently on the large oscillator. Therefore, a tran- 
sient QPO feature should appear around these times at a 
frequency close to ujq- In Fig. [8]and Fig. |9]we show the dy- 
namical spectrum from a simulation where the conditions 
for transient QPOs were engineered for two frequencies. 
The transient QPOs agree well with the expectations. As 
is seen from the figures, the strongest transient QPOs are 
those whose frequencies are the closest to that of the large 
oscillator; this is because the response of the large oscilla- 
tor is the strongest around its proper frequency. 

One can now easily understand the frequency drifts in 



Fig. 10 of L07 (Fig. 10 in this paper). In the simulations of 
that paper, the core continuum was sampled with a set of 
densely and regularly-placed Alfven modes by slicing the 
field into finite-width flux shells. The spacing Sco between 
the modes was not constant but a function of the Alfven 
frequency ui. In that case, the QPO drifts with the QPO 
frequency a;(i) given by the inverse relation 

With this relation we are able to fit all of L07 drifting 
QPOs, as shown in Fig. 10 and 11 Note that multiple 



QPOs correspond to different branches of the Alfven con- 
tinuum. As expected, the drifting QPOs amplified near the 
crustal frequencies, since there the response of the crust to 
the core modes' pull is the strongest. 
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Dynamical spectrum 




time 

Fig. 8. — Dynamical spectrum from a simulation where we 
have designed the continuum so as to produce transient QPO's 

at frequencies u = 1 and uj = 2 (the colored scale denotes Fig. 10. — Dynamical power spectrum of the spherical magne- 
log(power)). The green horizontal line denotes the frequency tar model from L07. The gray scale denotes log(power). 
of the large oscillator {Q — 1.2). 
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Fig. 9. — Dynamical spectrum from a simulation with a con- 
tinuum that is identical to the one of Fig ID We have shifted 
the frequency of the large oscillator, denoted by the green hor- 
izontal line, to n = 1.8. By comparison wit Fig. |8]it is clear 
that the drifting QPO's at oj = 2 are now much stronger as 
they are closer to the large oscillator frequency. Note that the 
edge mode at lj — 2.5 is clearly visible. 



Fig. 11.— We have used Eq. (|25| to fit the drifting QPO's 
from figure [lO] The red curves are n = 1 drifts, green curves 
are n = 2 and blue curves are n — 3. The higher frequency 
drifts originate from Alfven overtones. 
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4. More realistic magnetar models 

In this section we extend the constant magnetic field 
and constant-density magnetar model from LOT to include 
more realistic pressure and density profiles and more gen- 
eral (but still axisymmetric) magnetic field configurations. 
Our aim is to use this model to: 1) calculate numerically 
Alfven eigenmodes and their eigenfrequencies on different 
flux surfaces inside the star, in order to determine the 
continuous spectrum of the fluid core, and 2) use these 
modes to simulate the dynamics of a realistic magnetar. 
In order to calculate the Alfven eigenmodes and eigenfre- 
quencies for a realistic magnetar model, we employ the 
linearized equations of motion for an axisymmetric mag- 
netized, self-gravitating plasma. The general equations, 
which arc derived in detail in Poedts et al. (1985, hereafter 
P85) and given in their equations (53) and (54), consti- 
tute a fourth order system of coupled ordinary differential 
equations in the case of a mixed poloidal and toroidal mag- 
netic field. The formalism of P85 is briefly summarized 
in the Appendix B. In the case of a purely poloidal mag- 
netic field, the system simplifies to two uncoupled second 
order differential equations (P85, equations (70) and (71)). 



4.1. The model 

We assume our star is non-rotating and neglect its 
deformation due to the magnetic pressure, which is ex- 
pected to be small. Therefore, we consider a spheri- 
cally symmetric background model that is a solution of 
the Tolman-Oppcnheimcr-Volkoff equation (TOV equa- 
tion). The hydrostatic equilibrium is calculated using a 
SLy equation of state (Douchin & Haensel 2001; Haensel 
& Potekhin, 2004; Haensel, Potekhin & Yakovlev 2007), 
which can be found in tabulated form on the website 
hUp://www.ioffe.ru/astro/NSG/NSEOS/. The integra- 
tion of the TOV equation is performed using a 4th order 
Runge-Kutta scheme, integrating from the center of the 
star outward until we reach a mass density p = 1.3- 10^* 
g/cm^, which is consistent with the crust-core interface 
in the equation of state from Douchin & Haensel (2001). 
The resulting model has a central mass density pc = 10^^ 
g/cm^, a total mass of 1.40 Mq and a radius of i?coro = 
1.07- 10^ cm. To this spherical model we add a poloidal 
magnetic field, which we generate by placing a circular 
current loop of radius a and current I around the center 
of the star. The field is singular near the current loop, 
however all the field lines which connect to the crust (and 
thus are physically related to observable oscillations) carry 
finite field values. This particular field configuration is cho- 



sen as an example; there is an infinite number of ways to 
generate poloidal fleld configurations. In the appendix B 
we will add to this fleld a toroidal component and calculate 
the corresponding Alfven continuum of the core. 

4.2. The continuum 

In order to flnd the equations of motion for the magne- 
tized material in the neutron star core, we would need to 
add self-gravity to the ideal magnetohydrodynamic equa- 
tions. This problem has been solved by P85 in a tour the 
force mathematical approach. In that paper the authors 
assume a self-gravitating axisymmetric equilibrium with a 
field geometry consisting of mixed poloidal and toroidal 
field components, and they derive linearized equations of 
motion. For this fleld geometry it is convenient to work 
with so-called flux-coordinates {tjj,x,<l>)- The basic concept 
behind this curvilinear coordinate system is the magnetic 
flux-surface, which is deflned as the surface perpendicular 
to the Lorentz force Fl oc j x B. Prom this deflnition it 
is clear that the magnetic field lines lie in flux surfaces. If 
one considers a closed loop on a flux surface which makes 
one revolution around the axis of symmetry, then the mag- 
netic flux -0 through the loop depends on the flux surface 
only and is the same for all of the loops. Therefore ifi is 
chosen as the coordinate labeling the flux surfaces. In each 
flux-surface we can denote a position by its azimuthal an- 
gle (f) and its 'poloidal' coordinate X) which is deflned as 
the length along ^ = const line. In P85, it is shown that 
the equations of motion allow for a class of oscillatory so- 
lutions that are located on singular flux surfaces, consti- 
tuting a continuum of eigenmodes and eigenfrequencies. In 
the case of a purely poloidal field {B = B^), the continuum 
solutions are degenerate and polarized either parallel (^^) 
or perpendicular (^^) to the magnetic field lines. In the 
latter case the displacement is (^-independent. It is clear 
that in contrast to the x-polarized modes, the 0-polarized 
modes are purely horizontal and are therefore unaffected 
by gravity. This latter case is considered here. 

The equation of motion is in this case simply the Alfven 
wave equation: 



F[M,X)], 



where the operator F is given by 
B d 



FM^I^,X)\ = 



A-Kxp dx 



x'^B 



d ( Ui^,x) 

dx\ X 



(26) 



(27) 



Here x is the distance to the magnetic axis of symmetry. 
Although in the presence of a mixed poloidal and toroidal 
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field the equations still give rise to a continuous set of so- 
lutions, the calculations are significantly complicated as 
the continuum modes are affected by the toroidal compo- 
nent of the field, by gravity, and by compressibility. For 
the sake of simplicity we will ignore toroidal fields in our 
dynamic simulations. We will however, calculate the con- 
tinuum frequencies for a mixed poloidal and toroidal field 
in the Appendix B. 

For determining the spectrum of the core continuum, 
the appropriate boundary conditions arc ^^{x = Xc) = 0, 
where Xc{(t') marks the location of the crust-core inter- 
face. With this boundary condition. Equation (26) con- 



stitutes a Sturm-Liouville problem on each separate fiux 
surface ip. Using the stellar structure model and magnetic 
field configuration from section 4.1, we can calculate the 
eigenf unctions and eigenfrequencies for each flux surface 
tjj. The reflection symmetry of the stellar model and the 
magnetic field with respect to the equatorial plane assures 
that the eigenfunctions of Eq. ( 26 1 are either symmetric or 
anti-symmetric with respect to the equatorial plane. We 
can therefore determine the eigenfunctions by integrating 



Eq. ( 26 ) along the magnetic field lines from the equatorial 
plane x = to the crust-core interface x — Xc {4') ■ Let 
us consider the odd modes here for which (0) = 0, and 
solve Eq. (26) with the boundary condition ^0 (xc) = 
at the crust-core interface; for even modes, the boundary 
condition is d£^^ (0) /dx = 0. We find the eigenfunctions 
by means of a shooting method; using fourth order Runge- 
Kutta integration we integrate from X — ^ ^'^ X — Xc- The 
correct eigenvalues ct„ and eigenfunctions ^„ (%) are found 
by changing the value of a until the boundary condition at 
^„ is satisfied. In this way we gradually increase the value 
of a until the desired number of harmonics is obtained. In 
figure [12] we show a typical resulting core-continuum. 

According to Sturm-Liouville theory the normalized 
eigenfunctions ^„ of Eq. ( 26 1 form an orthonormal basis 



with respect to the following inner product: 



(Cm,Cn) = / r (x) £,m ix) ix) dx 



Where Sm n is the Kronecker delta and r 



(28) 



is the 



A7:p/B^ 

weight function. We have checked that the solutions we 
find satisfy the orthogonality relations. 

We are now ready to compute the coupled crust-core 
motion. Here we follow L07 and assume that the crust is 
an infinitely thin elastic shel]^ We label the lattitudinal 



^It is straightforward to relax this assumption, and carry out the anal- 
ysis of this section for the finite crustal thickness. However, from 
Section 2 it is clear that the interesting dynamics is dominated by 
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Fig. 12. — The red curves show the Alfven frequencies cr„ as 
a function of the angle ^('i/'), the polar angle at which the flux- 
surface 1/) intersects the crust. Since we are only considering 
odd crustal modes, the only Alfven modes that couple to the 
motion of the star are the ones with an odd harmonic number n. 
This particular continuum was calculated using a poloidal field 
with an average surface value Bsurfacc ~ 6- 10^'* G, generated 
by a circular ring current of radius a — R,/2. 
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Fig. 13. — After filling the curves from Fig. |12| 'gaps' in the 
continuum become visible around cr ~ 70 Hz and cr ~ 120 Hz. 
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location by the flux surface i/j intersecting the crust, and 
consider the crustal axisymmetric displacements ^^(V'')- In 
the MHD approximation, the magnetic stresses enforce a 
no-slip boundary condition at the crust-core interface, such 
that ^0 {ip, Xc) = £,4, {i', Xc) instead of (1/', Xc) = 0. It is 
useful to make the following substitution 



C (V", x) = ^0 (V", x) - (V') w (ip, x) 



(29) 



where we choose the function w x) so that (a) it cor- 
responds to the static displacement in the core and hence 
satisfies F{w{ip,x)) — 0, and (2) w (?/;, Xc) — 1- There- 
fore the new quantity satisfies the boundary condition 
C (V'j Xc) = and can be expanded into the Alfven normal 
modes ^„ which satisfy the same boundary conditions. 



We now proceed by substituting Eq. ( 29 1 into Eq. ( 26 ) 



thus obtaining a simple equation of motion for C 

75 ^(C(V',x)) = -w(^,x) — — (30) 



From the definition of the operator F it follows that for 
the odd modes 



u;(V',x) ^x{'il),x) 



xm,,x')B^{ij,x') 



r.dx'- (31) 



Here the constant K (ip) is chosen such that w {ip, Xc) = Ij 
in order that C = on both boundaries. We expand C and 
w into a series of ^„'s: 



\an {■>P,t)^n (V',X) 



(32) 

(33) 



Eq. (30) reduces to the following equations of motion for 



the eigenmode amplitudes a. 
d^an (ip) 



a<2 



'^n W O-n W = -bn W 



(34) 



These equations show how the core Alfven modes are 
driven by the motion of the crust. To close the system, 
we must address the motion of the crust driven by the 
hydromagnetic pull from the core. 

The equation of motion for the crust is given by 



92 



Lci (C0) + Lb 



(35) 



the spectral structure of the core Alfven waves; therefore in order to 
flesh out the physics we choose the simplified model of the crust. 



Where the acceleration due to elastic stresses Ld is 



COt(0)'-lU"0 



(36) 



where 6 is the polar angle (cf. L07). The acceleration Lb 
due to the magnetic stresses between the crust and the 
core can be expressed as 



xB^ d 

— cos a 7— 

47rS dx 



(37) 



where x is the distance to the axis of the star, S is column 
mass-density of the crust and a is the angle between the 
magnetic field line and the normal vector of the crust. 

It is convenient to express the crustal displacement £,4, as 
a Fourier series, being a sum normal modes of the free-crust 
problem. Using Eq. ( |36[ ) is straightforward to show analyt- 
ically that the eigenfunctions /; of the free-crust problem 



(Eq. (351 with Lb^O) are 



dYio {0) 



d9 



fi (0) oc 
with eigcnfrcquencies 



(38) 



(39) 



Here Yiq is the m — spherical harmonic of degree /. The 
normalized functions fi form an orthonormal basis, so that 



fi{0)U{e) sin{9)de = 6i. 



(40) 



where Si^m is again the Kronecker delta. The crustal dis- 
placement can then be expressed in terms of /; 



C4{9,t) = J2ci {t)fi{e) 



(41) 



Substituting Eq. (41 ) into Eq. (35 1 we obtain the equations 



of motion for the crustal mode amplitudes c/ 

Lb {9, t) fi (9) sin 9d9 



5^ 



(42) 



We can express Lb as 



Lb {ip,t) 



-^^cos(a(^)) 



a™ [t) a,. (43) 



dx 



K{^) 
X (7/;) B (^) 



X=Xc 
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Up to this point the derived equations of motion for the 
crust and the fluid core are exact. We are now ready to dis- 



cretize the continuum by converting the integral of Eq. ( 42 ) 
into a sum over N points 9i. In order to avoid the effect 
of phase coherence (see section 3) which caused drifts in 
the results from L07, we sample the continuum randomly 
over the 0-interval [0, 7r/2]. In the following, functional de- 
pendence of the coordinate %p as: 9 {i/j) is substituted by the 
discrete index i which denotes the i-th flux surface. 



L^f ci ^ 2 J2 L B i9r, t) fu sine, Ae, (44) 



= -^sin(^^,)A0./.i 



— cos(a,) L« 



Ckf ik 



dx 



X=Xc 



y air. 



I 



(45) 



These are the equations that fully describe dynamics of 
our magnetar model. As with the toy model from section 
2 we integrate them using a second order leap-frog scheme 
which conserves the total energy to high precision. As a 
test we keep track of the total energy of the system during 
the simulations. Further we have checked our results by 



integrating equations ( 45 1 and ( 45 ) with the fourth-order 



Runge-Kutta scheme and found good agreement with leap- 
frog integration. 

4.3. Results 

Based on our section-2 results, we have a good idea of 
what type of dynamical behavior should occur in our more 
realistic magnetar model. First, we expect that crustal 
modes with frequencies inside the Alfven continuum will 
be damped quickly by resonant absorption ("Landau- 
damping" in the terminology of Gruzinov 2008b). Second, 
as with our previous model we expect the late time be- 
havior of the system to show QPO's near the edges of the 
continuum, or edge modes. The third important feature 
of our model is that the continuum may possibly contain 
gaps, as is shown in Fig. [13] In this case there is the 
possibility that crustal frequencies fall inside the gaps and 
remain undamped. In all of our simulations these expecta- 
tions have come true. We will now show the results from 
a simulation which illustrate the above mentioned effects. 



The basic freedom of choice that we have for our model 
is the strength and geometry of the equilibrium magnetic 
field. We choose here a purely poloidal magnetic field with 
an average strength at the surface of Bsurf — 10^^ G, 
induced by a circular current loop of radius a = 0.5i?*. 
This field gives us a gap in the continuum at frequencies 
53 < w < 78 Hz. 

We consider the lowest degree odd crustal modes with 
frequencies W2 = 40 Hz and W4 — 84.5 Hz, which we cou- 
ple to 5000 continuum oscillators (the Alfven continuum). 
We sample the continuum at 1000 randomly chosen fiux 
surfaces, and at each flux surface we consider 5 Alfven 
overtones. 

Like with our toy model from section 2, we initiate the 
simulation by displacing the crust (c2 — — X) while 
keeping the continuum oscillators (the Alfven modes) re- 
laxed iain = 0). 

In Figures [14] and [15] we show the resulting power spec- 
trum for 2 different models. In the first one, the crustal 
frequencies are located inside the core continuum range, 
and the peaks due to the edge modes appear. By contrast, 
in the second case one of the crustal frequencies belongs 
to the gap, and a peak representing the global gap mode 
strongly stands our above the background. We note that 
the gap-mode's frequency lies close to but does not co- 
incide with the crustal-mode frequency; we found this to 
be a generic feature of our models, with the difference of 
10% for the typical model parameters. The gap modes are 
particularly interesting because they have relatively large 
amplitudes, and are not as strongly damped by viscosity 
as the edge modes. 

It must be emphasized that for all persistent modes in 
the system, the position in the frequency space of the core 
Alfven continuum plays the key role in setting the global- 
mode frequency and in determining its longevity. 

We note that Lee (2008) has used a different method to 
identify discrete modes in a magnetar with similar mag- 
netic configuration to ours. These modes were not asso- 
ciated with crustal frequencies, and we strongly suspect 
that they were located in the gaps of the continuum spec- 
trum and could be identified with the edge or gap modes 
presented in this work. 

5. Tangled magnetic fields 

Our preceding discussion of the continuum was predi- 
cated on the foliation of the axisymmetric magnetic field 
into the flux surfaces, with each of the singular contin- 
uum mode localized on the flux surfaces. These modes are 
"large" -they are coherent over the spacial extent compara- 
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Fig. 16. — Schematic illustration of tangled a magnetic field 
inside a magnetar. Locally, the field consists of flux tubes which 
contain a continuum of twisting Alfven modes. 



Fig. 14. — Power spectrum of the crustal dynamics for a mag- 
netar with a single 'gap' in the Alfven continuum. In this case 
the crustal frequencies are within the continuum, causing the 
crust modes to be Landau-damped. 
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Fig. 15. — Power spectrum of the crustal dynamics for a mag- 
netar with a single 'gap' in the Alfven continuum. The global 
mode within the gap is not damped, and its frequency is simi- 
lar, but not identical, to that of the crustal mode in the same 



ble to the size of the system, and thus they play an impor- 
tant role in the overall dynamics-they are responsible for 
the resonant absorption of the crust oscillations, and con- 
tribute to generating the edge and gap modes. But what 
happens if the field cannot be foliated into the flux sur- 
faces, but is instead tangled in a complicated way? One 
can argue that the continuum part of the spectrum still 
persists, as follows: 

Consider an arbitrary field line anchored at the crust- 
core interface at both ends, and choose a tube of field lines 
of infinitesimal radius which is centered on the original field 



line (see Fig. 16). It is clear that a twisting Alfven mode 



exists in the tube: it is obtained by the circular rotation 
of the fluid around the central field line, propagating along 
the central field line with the local Alfven velocity. Since 
there is a continuum of the field-line lengths, there is also 
a continuum of Alfven modes. However, the modes we 
constructed are highly localized in space and and have a 
small leverage in the overall dynamics. We conjecture that 
the more tangled the fields are, the less role do the singular 
continuum modes play in the overall dynamics. 

Whilst we cannot rigorously prove this conjecture, we 
can motivate it as follows: consider an area element SS of 
random orientation with the normal h inside the star, and 
consider a shearing motion along the element. This shear- 
ing motion will be resisted by the component of the 
magnetic field, with the effective shear modulus of order 



Mcff 



Itt 



(46) 



where (...) stands for averaging over the area element. For 
ordered field, it is possible to choose the orientation of the 
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area element so that /Zcs — 0; the presence of such an 
orientation makes a fundamental difference between MHD 
and elasticity theory and is responsible for the presence 
of continuous spectrum in MHD. However, if the linear 
size of the 5S is greater than the characteristic length on 
which the field is tangled, then /^eff is non-zero for all ori- 
entations of n. Therefore, for highly-tangled fields there 
can be no large-scale singular continuum modes, and their 
existence is restricted to the small scales. Hence our asser- 
tion that for strongly tangled fields continuum modes play 
a secondary dynamical role. 

One is then faced with the problem when crustal modes 
are coupled to a set of discrete core Alfven modes. In 
Appendix A we show how to find the eigen-solution of such 
a problem, provided that all of the coupling coefficients are 
known. 

How does one quantify the degree to which the fields 
are tangled? Some insight comes from the numerical sim- 
ulations of Braithwaite and colleagues, who have studied 
what type of fossil fields remain in a stratified star after 
an initial period of fast relaxation. Consider a stable fossil 
field field configuration, such as the one obtained in the 
simulations of Braithwaite and Spruit (2004) and Braith- 
waite and Nordlund (2006) [see also Gruzinov (2008a) for 
analytical considerations]. There, the final field is nearly, 
but not perfectly axisymmetric and has a small-scale ran- 
dom component. For a less-centrally concentrated initial 
field, Braithwaite (2008) shows that the final fossil field is 
in general non-axisymmetric and can have a complicated 
topology^ 

As a starting point, we shall consider the nearly axisym- 
metric field with a small random component. The latter 
acts like a small extra shear modulus ficB and dynamically 
couples the flux surfaces of the axisymmetric component. 
We then quantify the degree of tangling by the relative 
value of ficS and B'^/{4tt). 

5.1. simple model: "square" neutron star 

To study this idea further, we specify a very simple 
model of a neutron star, motivated by the one considered 



top plate 



in Levin (2006, hereafter L06) see Fig. 17 that never-the- 
less captures the essential physics. 

Consider a perfectly conducting homogeneous fluid of 
density p contained in a box with width , length Ly and 

^Gruzinov (2009) demonstrates that even this situation is not the most 
general. He finds that the relaxed field generally has multiple current 
sheets, and argues that the global field relaxation is dominated by the 
dissipation within these singular layers. The details do not concern 
us for the purposes of this paper. 
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Fig. 17. — Schematic illustration of the box model Per- 
fectly conducting incompressible fluid is sandwiched between 
perfectly conducting top and bottom plates. The box is peri- 
odic in 2-direction and the displacements of the plates (crust) 
are in the z-direction. The magnetic field is directed along y- 
axis and its strength varies as a function of x. 



depth Lz- The magnetic field in this box is everywhere 
aligned with the y-axis and its strength is a function of 
X only. We assume that gravity is zero and consider a 
Lagrangian displacement ^ (a;, y, t) of the fluid along the z- 
direction; we specify periodic boundary conditions in this 
direction (one should think of the z direction as azimuthal) . 
We now add to this model a small effective shear modulus 
/icff due to the field tangling. The fluid equation of motion 



9<2 



ca (a;) 



(47) 



Here ca (a;) is the Alfven velocity and Cg = /ieff/p is the 
/icff-generated shear velocity. If we assume a small shear 
speed, i.e. Cg « ca, Eq. (|47| reduces to 



(x) 



■dx^' 



(48) 



We now find the core Alfven eigenmodes. After adapting 
the no-slip boundary conditions 



0, 



(49) 



the problem can be easily solved by separation of vari- 
ables ^ {x, y,t) oc e"^' sin{7rm[(y/L^) + 1/2]} X (x), where 
TO — 1. 2, .... Plugging this in Eq. (48) we find for the the 
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a;— dependent part of the solution: 



2d'X _ 



(x)] X. 



(50) 



Here uJA.m i^) — Trmc^ (x) / Ly can be interpreted as the 
frequency of the m-th Alfven overtone at x. From the 
above expression it is clear that in the limit of very small 
Cs, the solution for X must be close to zero everywhere ex- 
cept in a very small neighborhood of ujA.m{x) = It is in 
this limit that the solutions are located on singular flux sur- 
faces. However, in the presence of the non- vanishing shear 
velocity Cg , the eigenmodes spread out on neighboring field 
lines, effectively coupling the motion on different flux sur- 
faces. The continuum of Alfven frequencies uiA,m ix) will 
in this case be no longer solutions of the system. Instead, 
the coupling term gives rise to a discrete set of solutions 
rather than a continuum. Eq. ( 50 ) is the mathematical 



equivalent of Schrodingers equation, which can in general 
cases be solved numerically. However, for many special 
case exact solutions exist. Let us consider, for the sake of 
simplicity, a field configuration in our box such that: 



c\ (x) 



2 I 2 
^X -t- 



We can rewrite Eq. (50) as follows: 



L,, 



-^x^X 



(51) 



X (52) 



This differential equation is the mathematical equivalent of 
the quantum harmonic oscillator problem for which the ex- 
act solution is well known. The eigenfrequencies are given 

by 



TT (1 + 2n) mCs 



'jLy + cloTT^m^/Ly 



(53) 



Here n {— 0,1,...) is the 'quantum' number labeling the 
harmonic-oscillator wavef unctions. We see that instead of 
a continuum, we obtain a densely packed discrete set of 
frequencies with the frequency spacing Wm.n ~ '^m,n-i ~ 

With the no-slip boundary conditions on the left and 
right sides x = ztLx/2, the eigenvalue equation must be 
solved numerically. An example of such calculation is 
shown in Fig. [l8] There, the spacing between the dis- 
crete Alfven modes is shown to increase as one increases 
the effect of the field tangling characterized by the fi^e- 

We now introduce the crustal modes into the problem by 
making the top and bottom of the box elastic and mobile. 
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Fig. 18. — Alfven frequencies as a function of the effective mag- 
netic shear modulus. As one decreases the shear, the spectrum 
tends to a continuum. 



We allow their displacement ^i.t,(a;,t) in the z-direction, 
and impose the boundary conditions on the sides: 



^tA-Lj2, t) = Ct,fc(i t) = 0. 



(54) 



Here the subscripts " i" and " 6" stand for the top and bot- 
tom of the box, respectively. The top and bottom are 
assumed to be thin and have mass M^r and surface den- 
sity fj = Met/ [LxLz) ■ The crustal equation of motion is 
given by 



9t2 



" dx^ Aira 

2^ , {BzBx}b 



dx'^ 



Ana 



(55) 



where Vs is the shear velocity in the crust. The crustal 
angular frequencies are given by — jirVg/Ly with 
the corresponding crustal mode eigenfunctions £,j = 
sin{jTT[{x/Lx) -1-1/2]} , where j — 1, 2, ... is roughly equiv- 
alent to I in the spherical case. The symmetry of the 
problem allows one to consider either symmetric — C& 
or antisymmetric = crustal modes. This will cou- 
ple to the symmetric (m = 1,3,5,...) or antisymmetric 
(to = 2, 4, 6, ...) Alfven modes of the core. 

Just as in section 4, it is now convenient to define a new 
variable C(x, y, t) for the core displacement: 



C{x, y, t) = ^{x, y, t) - £_o{x, y, t), 



(56) 
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where 



and 



y 



(57) 



The new variable observes the regular boundary condition 
C = on all the box edge, and satisfies the following inho- 
mogeneous partial differential equation: 



where 



9y2 



52 

^lg^]((.^^y^*)^9{x,y,t), (58) 



gix,y,t) 



52 \ 



(59) 



The advantage of the new variable is that it satisfies the 
regular boundary condition ^ = on all the boundaries of 
the box. It can therefore be expanded as a series consisting 
of eigenfunctions of the right-hand side of Equation 



(48): 



C(a;,y,i) = ran 

{x,y). (60) 

The rest of the procedure is very similar to that in sec- 
tion 4. We expand the crustal displacement into a series 
consisting of the eigenmode wavefunctions ^j: 

Ux,t) - J:,q,m,{x), (61) 

where Pj{t) and qj(t) are real numbers. The magnetar 
deformation is now fully represented by a set of generalized 
coordinates [p_,(t), gj(i), amn(i)]- The coupled equations of 
motion are derived by following the procedure specified in 
section 4. We obtain the following system of equations: 

2 



L,r, 



Pj 



1j 



(mn)j 



and 



Pj 



1j 



(63) 



^mn( 1) l^j{mn)^ran 



where 



a 



ip) 

(mn)j 
{mil)] 



/ (1 + L^) ^rnn{x,y)ij{x)dxdy 

mn 

/ (1 ^ L^) ^rnn{x,y)^j{x)dxdy 

mn 

{x,y)Ydxdy 



(64) 



dy 



V = Ly/2 



^j{x)dx 



(65) 



lUxWdx 

Thus we have obtained a system of linear second-order dif- 
ferential equations, which describes the time evolution of 
the square-box magnetar. These equations are solved by 
trancating all the series [i.e., rescricting the range of indices 
(to, n, j)] and then by either solving the eigenvalue prob- 
lem in order to find the normal modes, or by integrating 
the equations numericalljj^ One then checks that the se- 
ries trancation does not introduce errors in the magnetar's 
motion within the frequency range of our interest. 

So far we have worked in the approximation of the thin 
crust, i.e. we have effectively included the crustal modes 
which have no radial nodes in their wavefunction. How- 
ever, several observed high-frequency QPOs, and in par- 
ticular the strong QPO at 625Hz (Watts & Strohmayer 
2006) necessitate introduction of higher radial-order modes 
into our model. In the square-box model we do this phe- 
nomenologically, as follows. We assume that higher radial- 
order crustal modes have amplitudes Psj (t) and Qsj (t) and 
natural eigenfrequencies w^J, with s being the number of 
radial nodes, and assume that they cause displacement at 
the crust-core interface given by ^j{x). This mirrors real- 
istic spherically-symmetric case where the functional form 
of the crust-core displacement due to the torsional V x Y;,„ 
mode of the n'th radial order is a very weak function n. 
The amplitudes Psj{t) and qsj{t) are then introduced on 
into the equations of motions ( 62 1 and ( 63 1 in the same 



way as the other pj and qj amplitudes, with the same j- 
dependent coupling coefficients but with oj"' instead of 



on the left-hand side of Eq. ( 63 1 



We now have the basic ingredients of building a phe- 
nomenological modes with tangled fields. To sum up, we 
(1) quantify tangling using the effective shear modulus, (2) 
Find discrete core eigenmodes and evaluate their coupling 
to the crustal model, and (3) Either find eigenfrequencies 
of the total star by diagonalizing the potential energy of 
the system, or simulate the time-dependent behavior di- 
rectly. 

An example of a resulting powerspectrum is shown in 
Fig [19] for the model described in this section. 



'Our favored method here is again the energy-conserving second-order 
leapfrog. It is both fast and stable over long integration times. 
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Fig. 19. — Power spectrum for the dynamics of a magnetized 
box as described in the text. In this particular model we have 
used the maximum possible shear modulus, corresponding to 
a maximally tangled field. The Alfven motion in the box is 
coupled to 9 of the lowest frequency 'crustal' modes, plus a 
high frequency crust mode at 630 Hz. 



6. What do our models tell us about magnetar 
QPOs? 

In this paper we have developed a formaUsm which al- 
lows one to build a magnetar model with a variety of the 
spectral features of the core Alfven waves, including con- 
tinua with gaps and edges, and the large-scale discrete 
modes generated by the field tangling. We have con- 
structed a number of magnetar models and explored the 
resulting QPOS, both for the case of axisymmetric mag- 
netar with core Alfven continuum, and for the "square" 
magnetar models with the tangled fields (see the previous 
section). The full range of model parameters, and detailed 
comparison with the data will be the subject of a sepa- 
rate study. For now, we have restricted ourselves to the 
standard magnetar model, in which the core is a perfect 
conductor, the field of ^ 10^^ G penetrates both the core 
and the crust, and the proton fraction in the star is the 
one tabulated in Haensel, Potekin, and Yakovlev (2007). 
Our models give us the following robust conclusions, as 
compared against QPO observations: 

(1). A number of strong QPOs have been observed in 
the 1998 and 2004 giant flares, with frequencies ranging be- 
tween 18 Hz and 150 Hz ( Israel et al. 2005, Strohmayer and 
Watts 2006, Watts and Strohmayer 2006). These QPOs 
can be qualitatively explained as gap and/or edge modes 



Fig. 20. — This spectrum was generated using a box model 
similar to the one from figure[l9]but with neutron mass-loading. 
Due to the mass-loading the frequencies have shifted down by 
a factor of ~ 4. Note that there is no significant power above 
the lower edge-mode frequency of 5.3 Hz. 

of sections 4 and 2, or even transient QPOs of section ^ 
However, this was only possible if the neutrons were de- 
coupled from the Alfven waves in the core. If the neutrons 
took part in the Alfven motion, then the effective mass of 
the Alfven modes shifted up by a factor of 20 — 40 and 
their frequencies shifted down by a factor 4 — 8 (Easson 
& Pethick 1979, Alpar et al. 1984, van Hoven & Levin 
2008, Andersson et al. 2009). As a result, all modes at 
frequencies above ^ 50Hz were strongly damped (see Fig. 
20). Increasing the magnetic- field tension by a factor of 3 
did not affect this conclusion (Fig 21 1. For the spherical 
magnetar models of section 4 we obtain similar results if 
couple the neutrons to the Alfven motion in the core. The 
key point that we would like the reader to appreciate is 
that Alfven modes in the core are key to determining the 
frequency and strength of the observable QPOs, and thus 
QPOs are very sensitive probe of the core interior. 

(2) A number of the high-frequency QPOs have been 
measured in the 2004 giant flare by Watts and Strohmayer 



^L07 and Gruzinov 2008b associated the long-lived 18-20Hz QPO with 
the lower edge of the Alfven continuum. However, recent calculations 
of Steiner and Watts (2009) have argued that the crustal frequen- 
cies can be as low as lOHz due to the uncertainty in our theoretical 
knowledge of the crustal shear modulus. It is therefore plausible that 
the fundamental crustal mode has the proper frequency below that 
lower edge of the core Alfven continuum. In this case, the 18-20Hz 
QPO could be the gap mode which is dominated by the fundamental 
crustal mode. 
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Fig. 21. — This spectrum was generated with the same box 
model as in figures |20[ but in addition to the neutron mass- 
loading, we have increased the magnetic field strength by a 
factor of 3. All frequencies above ~ 16 Hz are significantly 
damped. 



(2006), the strongest among them being the QPO at 625 
Hz. This QPO is particularly strong and long-lived in the 
hard x-rays, reaching the amplitude of ^ 25% over the time 
interval of ~ 100 seconds (i.e., it persists for almost 10^ os- 
cillation periods!). Watts and Strohmayer (2006) argued 
that this frequency corresponds to the crustal shear mode 
with a single radial node (see also Piro 2005) ; this interpre- 
tation, if correct, would strongly constrain the thickness of 
the crust and rule out the fluid strange stars as magne- 
tar candidates (Watts & Reddy, 2007). To investigate this 
suggestion, we have introduced several high-frequency low- 
j crustal modes into our square-box simulations. However, 
as is demonstrated in Figs. [T4]and[T5j the high-frequency 
modes are strongly damped and at no time during the 
simulations do we observe any significant power at those 
frequencies. This is to be expected. No natural axisym- 
metric model has gaps in the Alfven continuum at such 
high frequencies, so global modes are strongly absorbed. 
One could expect that if the Alfven modes are discrete 
in the core due to the field tangling, this problem would 
not arise. However, even in the discrete case the frequency 
spacing between the modes is around 20 Hz, which is much 
smaller than 600 Hz. Thus the grid of Alfven waves is so 
dense that it is effectively seen as the absorbing continuum 
by the modes around 600 Hz. Our detailed simulations, of 



the type shown in Figs. 14 and 15 fully confirm this qual- 
itative picture. 



The concern about the viability of high-frequency QPOs 
as being due to the physical oscillations of standard-model 
magnetars has been raised in the original L06 paper on the 
basis of rather simplistic calculations. As our work here 
shows, more detailed calculations partially alleviate the 
L06 concern for the frequencies below ~ 150Hz, but only 
if the neutrons are decoupled from the Alfven motion in 
the core, i.e. if at least one baryonic superfluid (protons or 
neutrons) are present in the neutron-star core. Our anal- 
ysis sustains L06 concern for the high-frequency QPOs, in 
particular for the strong long-lived QPO at 625Hz. Its ex- 
planation seems to require either QPO production in the 
magnetosphere, or a somewhat radical revision of the mag- 
netar model. Just how radical this revision has to be will 
be explored in a separate study. 

Our work presented here has several shortcomings. We 
have limited ourselves to the linear approximation, and a 
non-linear regime may bring surprises. Direct non-linear 
simulations of axisymmetric oscillations of a magnetised 
fluid star has been carried out recently by Cerda-Duran, 
Stergioulas, & Font (2009). At this stage it is difficult 
to say whether non-linearities introduce significantly new 
QPO features to their model; their results have largely 
been in agreement with the linear simulations of Colaiuda, 
Beyer, & Kokkotas (2009). However, the computational 
techniques seem promising and we do not exclude that 
large-amplitude simulations of stars with the crust will 
show qualitatively new features. Another limitation of our 
work is that we have assumed that once the flare sets the 
magnetar into motion, the magnetar's oscillations are not 
driven externally. This may not be the case in real flares: 
some energy stored in the pre-flare magnetar may be re- 
leased gradually, and this release could be extended in time 
into the flare's taij^ The latter consideration is straight- 
forward to incorporate phenomenologically into our model, 
and we plan to address it in our future work. 
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A. Multimodal crust-core system 

In this Appendix we generalize the normal-mode treatment of Section 2.2, and write down the general prescription 
of how to find the eigenmodes when several "large" crustal shear modes are coupled to a multitude of small core Alfven 
modes, provided the coupling coefficients are known. In this paper, the coupling coefficients are worked out in simple 
models of sections 4 and 5; we shall postpone the discussion of how the coefficients are computed in more general case 
to the future paper. 

Let us denote the displacement of the crustal and core modes by A„ and Xi respectively. Since both the crustal 
and the core modes are not directly coupled to themselves (i.e., X's are only coupled to a;'s), most general equations of 
motion take the form 



Xn + ri^A„ = T^iUniX^ (Al) 

where r2„ and uij are the proper frequencies of the crustal and core modes, and a's and /?'s are the coupling coefficients. 
We look for an oscillatory solutions of the above equations with angular frequency fl. One can trivially re-write these 
equation as a matrix eigen-equation with as an eigenvalue, and solve it using standard methods. However, if the 



number of crustal modes is not too large, it is convenient to make a shortcut. Using the second of Eq. ( Al ) to express 
Xi's through X„'s, and substituting into the first one, we get the following equation: 

S„G„„,(f])A„ = 0, (A2) 

where the elements of the matrix G are given by 

G™„(f7) = (n^ - nl)s„,r^ + s,^^. (A3) 

One obtains the eigenfrequencies by finding numerically the zeros of det Gmn- 



B. Core continua with a mixed axisymmetric toroidal-poloidal magnetic field 

In this appendix we will calculate the continuum of Alfven frequencies in a magnetar core in the case of a axisymmetric 
magnetic field with mixed toroidal and poloidal components. The general MHD equations of motion for spherically 
symmetric, self-gravitating equilibrium with an axisymmetric field, are derived in detail in P85. In contrast to the special 
case of a purely poloidal field (see section 4.2) which leads to two uncoupled differential equations, the continuum for a 
mixed toroidal-poloidal field is described by a system of fourth order coupled ODEs. Due to this coupling, the solutions 
are complicated as they are no longer polarized in the directions parallel (so-called "cusp solutions") and perpendicular 
(Alfven solutions) to the magnetic field lines, but rather have a mixed character. Strictly speaking, one can only speak 
of an "Alfven continuum" in the limit that the variations in p, P and are small in the magnetic flux-surfaces. The 
general equations of motion are given in Eqs. (53) and (54) of P85. We note however, that in magnetars the speed of 
sound c >> CA, and therefore we consider Poedts et al.'s equations (53) and (54) in the incompressible limit (P85, Eqs. 
(73) and (74)). For completeness we give the equations here, 

Y + pBlNl{Y + Z)- — {pc\)FZ (Bl) 

+ pBlNl {Y + Z) + F {pc\FZ) (B2) 

The variables Y = i i^B'^S.x ^ Bij^B^^^j /B^B^ and Z = i (B^^x + -^0^0) /-^^ components of the fluid displace- 
ment perpendicular and parallel to the magnetic field lines, the operator F = id/dx is a differential operator along 



pa 



2B^y 
Bl 



B^F 



{pel) 



pa^B'^Z 



iF 



dx 
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Fig. 22. — The curves show the continuum frequencies ct„ as a function of the angle the polar angle at which the flux-surface 

tp intersects the crust. In the presence of a toroidal field, the degeneracy between the cusp-solutions and the Alfven solutions is 
broken and we find two separate solutions for each wave number n; waves with primarily Alfven character (red curves) and waves 
with primarily cusp character (blue curves). This particular continuum was calculated using a poloidal field with an average 
surface value Bp, surface ~ 6- 10^* G (again generated by a circular ring curr ent of radius a = -R*/2) and a toroidal field strength 



at the equator and the crust-core interface of -Bt.o = 3-10 G (see Eq. (B3l) 



the field lines, = - {l/B^p) y/{dp/dx) (dP/dx) can be thought of as a Brunt- Vaisala frequency for displacements 
along the field lines. According to Gauss' law for magnetism, the toroidal component of the magnetic field is of the 
form = f{ip)/w^ where tu is the distance to the polar axis, and f{ip) is an arbitrary function of tjj. In the following 
calculation we adopt a toroidal field component of the form 

Here 9{ip) is the polar angle at which the flux surface ip intersects the stellar crust. Clearly this choice for i?^ is 
completely arbitrary and one could in principle try many different toroidal geometries. 

As with our calculation of the Alfven continuum in the case of a purely poloidal field (section 4.2), we adopt the zero- 
displacement boundary conditions at the crust, and use the fact that our equilibrium model is (point-) symmetric with 
respect to the equatorial plane. This enforces the existence of classes of symmetric and anti-symmetric eigenfunctions, 
Ynix) and Zn{x)- We consider only the odd modes, and calculate the eigenfunctions by means of the shooting method; 



we use a fourth order Runge-Kutta scheme to integrate Eqs. (Bl) and (B2). Starting with Y{0) = and Z{0) = at 
the equator, we integrate outward until we reach the crust Bt x — Xc- We find the eigenfrequencies by changing the 
value of fj until we match the boundary conditions at the crust. A resulting continuum is plotted in Figure [22] 
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